Contents

1 Load dependencies

library(SingleCellExperiment)
library(BiocParallel)
library(geneBasisR)
library(splatter)
library(scran)
library(miloR)
library(miloDE)
library(stringr)
library(ggplot2)
library(ggpubr)
library(viridis)
library(MetBrewer)
library(dplyr)
library(tibble)

ncores = 24
mcparam = MulticoreParam(workers = ncores)
register(mcparam)
set.seed(32)

root.dir = "/nfs/research/marioni/alsu/hubmap_metaRef/"
source(paste0(root.dir , "miloDE_analysis/core_functions.R"))
figures.dir = paste0(root.dir , "figures/da_de/simulations_splatter/random/")

# generated simulation data
sces = readRDS(file = paste0(root.dir , "data/processed/da_de/simulations_random/sims.Rds"))
simulated_genes = readRDS(file = paste0(root.dir , "data/processed/da_de/simulations_random/simulated_genes.Rds"))
stat_da = readRDS(file = paste0(root.dir , "data/processed/da_de/simulations_random/stat_da.Rds"))
stat_de = readRDS(file = paste0(root.dir , "data/processed/da_de/simulations_random/stat_de.Rds"))

stat_combined = readRDS(file = paste0(root.dir , "data/processed/da_de/simulations_random/stat_combined_all_genes.Rds"))

2 Heatmaps

2.1 DA

2.1.1 Split by all

stat_combined_reduced = unique(stat_combined[, c("id_sim" , "id_nhood_assignment" , "reducedDim_name", "frac_hoods_DA_001", "frac_hoods_DA_005", "frac_hoods_DA_01" , "n_markers" , "n_not_markers" , "round_sim" , "order" , "k" , "seed")] )
stat_combined_reduced$order_k = paste0(stat_combined_reduced$order , "_" , stat_combined_reduced$k)

stat_combined_reduced$reducedDim_name[stat_combined_reduced$reducedDim_name == "pca.hvgs"] = "Supervised"
stat_combined_reduced$reducedDim_name[stat_combined_reduced$reducedDim_name == "pca.all"] = "Unsupervised"

get_plot = function(var){
  p = ggplot(stat_combined_reduced , aes(x = factor(n_markers) , y = factor(n_not_markers) , fill = stat_combined_reduced[,var])) +
    geom_tile() + 
    facet_wrap(~reducedDim_name + round_sim + order_k , ncol = 12) + 
    scale_fill_viridis(discrete = F , name = "Fraction of DA nhoods") +
    theme_bw() + 
    theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) + 
    labs(x = "# of altered genes contributing to joint variance" , y = "# of altered genes contributing to case variance") +
    theme(legend.position = "top")
  p
  return(p)
}


p = get_plot("frac_hoods_DA_01")
ggsave(filename = paste0(figures.dir, "da/split_all__fdr_01", ".png"), plot = p, width = 20, height = 6)

p = get_plot("frac_hoods_DA_005")
ggsave(filename = paste0(figures.dir, "da/split_all__fdr_005", ".png"), plot = p, width = 20, height = 6)

p = get_plot("frac_hoods_DA_001")
ggsave(filename = paste0(figures.dir, "da/split_all__fdr_001", ".png"), plot = p, width = 20, height = 6)

2.1.2 Groupped by same sim

stat_combined_reduced = unique(stat_combined[, c("id_sim" , "id_nhood_assignment" , "reducedDim_name", "frac_hoods_DA_001", "frac_hoods_DA_005", "frac_hoods_DA_01" , "n_markers" , "n_not_markers" , "round_sim" , "order" , "k" , "seed")] )
stat_combined_reduced$order_k = paste0(stat_combined_reduced$order , "_" , stat_combined_reduced$k)

stat_combined_reduced$reducedDim_name[stat_combined_reduced$reducedDim_name == "pca.hvgs"] = "Supervised"
stat_combined_reduced$reducedDim_name[stat_combined_reduced$reducedDim_name == "pca.all"] = "Unsupervised"

stat_combined_reduced = as.data.frame(stat_combined_reduced %>% group_by(id_sim , reducedDim_name, frac_hoods_DA_001, frac_hoods_DA_005, frac_hoods_DA_01 , n_markers , n_not_markers , round_sim) %>% dplyr::summarise(frac_hoods_DA_001 = mean(frac_hoods_DA_001 , na.rm = T) , frac_hoods_DA_005 = mean(frac_hoods_DA_005 , na.rm = T) , frac_hoods_DA_01 = mean(frac_hoods_DA_01 , na.rm = T)) )



get_plot = function(var){
  p = ggplot(stat_combined_reduced , aes(x = factor(n_markers) , y = factor(n_not_markers) , fill = stat_combined_reduced[,var])) +
    geom_tile() + 
    facet_wrap(~reducedDim_name + round_sim) + 
    scale_fill_viridis(discrete = F , name = "Fraction of DA nhoods") +
    theme_bw() + 
    theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) + 
    labs(x = "# of altered genes contributing to joint variance" , y = "# of altered genes contributing to case variance") +
    theme(legend.position = "top")
  p
  return(p)
}


p = get_plot("frac_hoods_DA_01")
ggsave(filename = paste0(figures.dir, "da/groupped_sim__fdr_01", ".png"), plot = p, width = 8, height = 6)

p = get_plot("frac_hoods_DA_005")
ggsave(filename = paste0(figures.dir, "da/groupped_sim__fdr_005", ".png"), plot = p, width = 8, height = 6)

p = get_plot("frac_hoods_DA_001")
ggsave(filename = paste0(figures.dir, "da/groupped_sim__fdr_001", ".png"), plot = p, width = 8, height = 6)

2.1.3 Groupped by eveyrhting

stat_combined_reduced = unique(stat_combined[, c("id_sim" , "id_nhood_assignment" , "reducedDim_name", "frac_hoods_DA_001", "frac_hoods_DA_005", "frac_hoods_DA_01" , "n_markers" , "n_not_markers" , "round_sim" , "order" , "k" , "seed")] )
stat_combined_reduced$order_k = paste0(stat_combined_reduced$order , "_" , stat_combined_reduced$k)

stat_combined_reduced$reducedDim_name[stat_combined_reduced$reducedDim_name == "pca.hvgs"] = "Supervised"
stat_combined_reduced$reducedDim_name[stat_combined_reduced$reducedDim_name == "pca.all"] = "Unsupervised"

stat_combined_reduced = as.data.frame(stat_combined_reduced %>% group_by(reducedDim_name, frac_hoods_DA_001, frac_hoods_DA_005, frac_hoods_DA_01 , n_markers , n_not_markers ) %>% dplyr::summarise(frac_hoods_DA_001 = mean(frac_hoods_DA_001 , na.rm = T) , frac_hoods_DA_005 = mean(frac_hoods_DA_005 , na.rm = T) , frac_hoods_DA_01 = mean(frac_hoods_DA_01 , na.rm = T)) )



get_plot = function(var){
  p = ggplot(stat_combined_reduced , aes(x = factor(n_markers) , y = factor(n_not_markers) , fill = stat_combined_reduced[,var])) +
    geom_tile() + 
    facet_wrap(~reducedDim_name) + 
    scale_fill_viridis(discrete = F , name = "Fraction of DA nhoods") +
    theme_bw() + 
    theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) + 
    labs(x = "# of altered genes contributing to joint variance" , y = "# of altered genes\ncontributing to case variance") +
    theme(legend.position = "top")
  p
  return(p)
}


p = get_plot("frac_hoods_DA_01")
ggsave(filename = paste0(figures.dir, "da/groupped_all__fdr_01", ".png"), plot = p, width = 5, height = 3)

p = get_plot("frac_hoods_DA_005")
ggsave(filename = paste0(figures.dir, "da/groupped_all__fdr_005", ".png"), plot = p, width = 5, height = 3)

p = get_plot("frac_hoods_DA_001")
ggsave(filename = paste0(figures.dir, "da/groupped_all__fdr_001", ".png"), plot = p, width = 5, height = 3)

2.2 DE - overall

2.2.1 Split by all

stat_combined_reduced = unique(stat_combined[, c("id_sim" , "id_nhood_assignment" , "reducedDim_name", "gene_type", "frac_hoods_DE_001", "frac_hoods_DE_005", "frac_hoods_DE_01" , "n_markers" , "n_not_markers" , "round_sim" , "order" , "k" , "seed")] )
stat_combined_reduced$order_k = paste0(stat_combined_reduced$order , "_" , stat_combined_reduced$k)

stat_combined_reduced$reducedDim_name[stat_combined_reduced$reducedDim_name == "pca.hvgs"] = "Supervised"
stat_combined_reduced$reducedDim_name[stat_combined_reduced$reducedDim_name == "pca.all"] = "Unsupervised"

get_plot = function(var , gene_type){
  current.stat = stat_combined_reduced[stat_combined_reduced$gene_type == gene_type , ]
  p = ggplot(current.stat , aes(x = factor(n_markers) , y = factor(n_not_markers) , fill = current.stat[,var])) +
    geom_tile() + 
    facet_wrap(~reducedDim_name + round_sim + order_k , ncol = 12) + 
    scale_fill_viridis(discrete = F , name = "Fraction of DE nhoods") +
    theme_bw() + 
    theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) + 
    labs(x = "# of altered genes contributing to joint variance" , y = "# of altered genes contributing to case variance") +
    theme(legend.position = "top")
  p
  return(p)
}


anno = expand.grid(gene_type = c("marker" , "not_marker") , var = c("01" , "005" , "001"))


plots = lapply(1:nrow(anno) , function(i){
  p = get_plot(paste0("frac_hoods_DE_" ,as.character(anno$var[i])) , as.character(anno$gene_type[i]))
  ggsave(filename = paste0(figures.dir, "de/" , as.character(anno$gene_type[i]), "/split_all__fdr_" , as.character(anno$var[i]), ".png"), plot = p, width = 20, height = 6)
})

2.2.2 Groupped by same sim

stat_combined_reduced = unique(stat_combined[, c("id_sim" , "id_nhood_assignment" , "reducedDim_name", "gene_type", "frac_hoods_DE_001", "frac_hoods_DE_005", "frac_hoods_DE_01" , "n_markers" , "n_not_markers" , "round_sim" , "order" , "k" , "seed")] )
stat_combined_reduced$order_k = paste0(stat_combined_reduced$order , "_" , stat_combined_reduced$k)

stat_combined_reduced$reducedDim_name[stat_combined_reduced$reducedDim_name == "pca.hvgs"] = "Supervised"
stat_combined_reduced$reducedDim_name[stat_combined_reduced$reducedDim_name == "pca.all"] = "Unsupervised"


stat_combined_reduced = as.data.frame(stat_combined_reduced %>% group_by(id_sim , reducedDim_name, gene_type , n_markers , n_not_markers , round_sim) %>% dplyr::summarise(frac_hoods_DE_001 = mean(frac_hoods_DE_001 , na.rm = T) , frac_hoods_DE_005 = mean(frac_hoods_DE_005 , na.rm = T) , frac_hoods_DE_01 = mean(frac_hoods_DE_01 , na.rm = T)) )



get_plot = function(var , gene_type){
  current.stat = stat_combined_reduced[stat_combined_reduced$gene_type == gene_type , ]
  p = ggplot(current.stat , aes(x = factor(n_markers) , y = factor(n_not_markers) , fill = current.stat[,var])) +
    geom_tile() + 
    facet_wrap(~reducedDim_name + round_sim ) + 
    scale_fill_viridis(discrete = F , name = "Fraction of DE nhoods") +
    theme_bw() + 
    theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) + 
    labs(x = "# of altered genes contributing to joint variance" , y = "# of altered genes contributing to case variance") +
    theme(legend.position = "top")
  p
  return(p)
}


anno = expand.grid(gene_type = c("marker" , "not_marker") , var = c("01" , "005" , "001"))


plots = lapply(1:nrow(anno) , function(i){
  p = get_plot(paste0("frac_hoods_DE_" ,as.character(anno$var[i])) , as.character(anno$gene_type[i]))
  ggsave(filename = paste0(figures.dir, "de/" , as.character(anno$gene_type[i]), "/groupped_sim__fdr_" , as.character(anno$var[i]), ".png"), plot = p, width = 8, height = 6)
})

2.2.3 Groupped by eveyrhting

stat_combined_reduced = unique(stat_combined[, c("id_sim" , "id_nhood_assignment" , "reducedDim_name", "gene_type", "frac_hoods_DE_001", "frac_hoods_DE_005", "frac_hoods_DE_01" , "n_markers" , "n_not_markers" , "round_sim" , "order" , "k" , "seed")] )
stat_combined_reduced$order_k = paste0(stat_combined_reduced$order , "_" , stat_combined_reduced$k)

stat_combined_reduced$reducedDim_name[stat_combined_reduced$reducedDim_name == "pca.hvgs"] = "Supervised"
stat_combined_reduced$reducedDim_name[stat_combined_reduced$reducedDim_name == "pca.all"] = "Unsupervised"


stat_combined_reduced = as.data.frame(stat_combined_reduced %>% group_by(reducedDim_name, gene_type , n_markers , n_not_markers , round_sim) %>% dplyr::summarise(frac_hoods_DE_001 = mean(frac_hoods_DE_001 , na.rm = T) , frac_hoods_DE_005 = mean(frac_hoods_DE_005 , na.rm = T) , frac_hoods_DE_01 = mean(frac_hoods_DE_01 , na.rm = T)) )



get_plot = function(var , gene_type){
  current.stat = stat_combined_reduced[stat_combined_reduced$gene_type == gene_type , ]
  p = ggplot(current.stat , aes(x = factor(n_markers) , y = factor(n_not_markers) , fill = current.stat[,var])) +
    geom_tile() + 
    facet_wrap(~reducedDim_name ) + 
    scale_fill_viridis(discrete = F , name = "Fraction of DE nhoods" ) +
    theme_bw() + 
    theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) + 
    labs(x = "# of altered genes contributing to joint variance" , y = "# of altered genes\ncontributing to case variance") +
    theme(legend.position = "top")
  p
  return(p)
}


anno = expand.grid(gene_type = c("marker" , "not_marker") , var = c("01" , "005" , "001"))


plots = lapply(1:nrow(anno) , function(i){
  p = get_plot(paste0("frac_hoods_DE_" ,as.character(anno$var[i])) , as.character(anno$gene_type[i]))
  ggsave(filename = paste0(figures.dir, "de/" , as.character(anno$gene_type[i]), "/groupped_all__fdr_" , as.character(anno$var[i]), ".png"), plot = p, width = 5, height = 3)
})

2.3 DE within not DA

2.3.1 Split by all

stat_combined_reduced = unique(stat_combined[, c("id_sim" , "id_nhood_assignment" , "reducedDim_name", "gene_type", "frac_hoods_DE_not_DA_001", "frac_hoods_DE_not_DA_005", "frac_hoods_DE_not_DA_01" , "n_markers" , "n_not_markers" , "round_sim" , "order" , "k" , "seed")] )
stat_combined_reduced$order_k = paste0(stat_combined_reduced$order , "_" , stat_combined_reduced$k)

stat_combined_reduced$reducedDim_name[stat_combined_reduced$reducedDim_name == "pca.hvgs"] = "Supervised"
stat_combined_reduced$reducedDim_name[stat_combined_reduced$reducedDim_name == "pca.all"] = "Unsupervised"

get_plot = function(var , gene_type){
  current.stat = stat_combined_reduced[stat_combined_reduced$gene_type == gene_type , ]
  p = ggplot(current.stat , aes(x = factor(n_markers) , y = factor(n_not_markers) , fill = current.stat[,var])) +
    geom_tile() + 
    facet_wrap(~reducedDim_name + round_sim + order_k , ncol = 12) + 
    scale_fill_viridis(discrete = F , name = "Fraction of DE nhoods\n(in not DA nhoods)") +
    theme_bw() + 
    theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) + 
    labs(x = "# of altered genes contributing to joint variance" , y = "# of altered genes contributing to case variance") +
    theme(legend.position = "top")
  p
  return(p)
}


anno = expand.grid(gene_type = c("marker" , "not_marker") , var = c("01" , "005" , "001"))


plots = lapply(1:nrow(anno) , function(i){
  p = get_plot(paste0("frac_hoods_DE_not_DA_" ,as.character(anno$var[i])) , as.character(anno$gene_type[i]))
  ggsave(filename = paste0(figures.dir, "de_in_not_da/" , as.character(anno$gene_type[i]), "/split_all__fdr_" , as.character(anno$var[i]), ".png"), plot = p, width = 20, height = 6)
})

2.3.2 Groupped by same sim

stat_combined_reduced = unique(stat_combined[, c("id_sim" , "id_nhood_assignment" , "reducedDim_name", "gene_type", "frac_hoods_DE_not_DA_001", "frac_hoods_DE_not_DA_005", "frac_hoods_DE_not_DA_01" , "n_markers" , "n_not_markers" , "round_sim" , "order" , "k" , "seed")] )
stat_combined_reduced$order_k = paste0(stat_combined_reduced$order , "_" , stat_combined_reduced$k)

stat_combined_reduced$reducedDim_name[stat_combined_reduced$reducedDim_name == "pca.hvgs"] = "Supervised"
stat_combined_reduced$reducedDim_name[stat_combined_reduced$reducedDim_name == "pca.all"] = "Unsupervised"


stat_combined_reduced = as.data.frame(stat_combined_reduced %>% group_by(id_sim , reducedDim_name, gene_type , n_markers , n_not_markers , round_sim) %>% dplyr::summarise(frac_hoods_DE_001 = mean(frac_hoods_DE_not_DA_001 , na.rm = T) , frac_hoods_DE_005 = mean(frac_hoods_DE_not_DA_005 , na.rm = T) , frac_hoods_DE_01 = mean(frac_hoods_DE_not_DA_01 , na.rm = T)) )



get_plot = function(var , gene_type){
  current.stat = stat_combined_reduced[stat_combined_reduced$gene_type == gene_type , ]
  p = ggplot(current.stat , aes(x = factor(n_markers) , y = factor(n_not_markers) , fill = current.stat[,var])) +
    geom_tile() + 
    facet_wrap(~reducedDim_name + round_sim ) + 
    scale_fill_viridis(discrete = F , name = "Fraction of DE nhoods\n(in not DA nhoods)") +
    theme_bw() + 
    theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) + 
    labs(x = "# of altered genes contributing to joint variance" , y = "# of altered genes contributing to case variance") +
    theme(legend.position = "top")
  p
  return(p)
}


anno = expand.grid(gene_type = c("marker" , "not_marker") , var = c("01" , "005" , "001"))


plots = lapply(1:nrow(anno) , function(i){
  p = get_plot(paste0("frac_hoods_DE_" ,as.character(anno$var[i])) , as.character(anno$gene_type[i]))
  ggsave(filename = paste0(figures.dir, "de_in_not_da/" , as.character(anno$gene_type[i]), "/groupped_sim__fdr_" , as.character(anno$var[i]), ".png"), plot = p, width = 8, height = 6)
})

2.3.3 Groupped by eveyrhting

stat_combined_reduced = unique(stat_combined[, c("id_sim" , "id_nhood_assignment" , "reducedDim_name", "gene_type", "frac_hoods_DE_not_DA_001", "frac_hoods_DE_not_DA_005", "frac_hoods_DE_not_DA_01" , "n_markers" , "n_not_markers" , "round_sim" , "order" , "k" , "seed")] )
stat_combined_reduced$order_k = paste0(stat_combined_reduced$order , "_" , stat_combined_reduced$k)

stat_combined_reduced$reducedDim_name[stat_combined_reduced$reducedDim_name == "pca.hvgs"] = "Supervised"
stat_combined_reduced$reducedDim_name[stat_combined_reduced$reducedDim_name == "pca.all"] = "Unsupervised"


stat_combined_reduced = as.data.frame(stat_combined_reduced %>% group_by(reducedDim_name, gene_type , n_markers , n_not_markers , round_sim) %>% dplyr::summarise(frac_hoods_DE_001 = mean(frac_hoods_DE_not_DA_001 , na.rm = T) , frac_hoods_DE_005 = mean(frac_hoods_DE_not_DA_005 , na.rm = T) , frac_hoods_DE_01 = mean(frac_hoods_DE_not_DA_01 , na.rm = T)) )



get_plot = function(var , gene_type){
  current.stat = stat_combined_reduced[stat_combined_reduced$gene_type == gene_type , ]
  p = ggplot(current.stat , aes(x = factor(n_markers) , y = factor(n_not_markers) , fill = current.stat[,var])) +
    geom_tile() + 
    facet_wrap(~reducedDim_name ) + 
    scale_fill_viridis(discrete = F , name = "Fraction of DE nhoods" ) +
    theme_bw() + 
    theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) + 
    labs(x = "# of altered genes contributing to joint variance" , y = "# of altered genes\ncontributing to case variance") +
    theme(legend.position = "top")
  p
  return(p)
}


anno = expand.grid(gene_type = c("marker" , "not_marker") , var = c("01" , "005" , "001"))


plots = lapply(1:nrow(anno) , function(i){
  p = get_plot(paste0("frac_hoods_DE_" ,as.character(anno$var[i])) , as.character(anno$gene_type[i]))
  ggsave(filename = paste0(figures.dir, "de_in_not_da/" , as.character(anno$gene_type[i]), "/groupped_all__fdr_" , as.character(anno$var[i]), ".png"), plot = p, width = 5, height = 3)
})

2.4 DA vs DE

@@ Mb later, lets see if we need it

3 UMAP cartoons

3.1 Plotting function

cols = met.brewer("Egypt",4)
cols = cols[c(3:4)]

get_plot = function(id , red_dim){
  sce = sces[[id]]

  umaps = as.data.frame( scater::calculateUMAP(t(reducedDim(sce , red_dim))) )
  reducedDim(sce , "umap") = umaps
  meta = cbind(as.data.frame(colData(sce)) , as.data.frame(reducedDim(sce , "umap")))
  meta$Group = as.character(meta$Group)
  meta$Group[meta$Group == "Group1"] = "Cell type 1"
  meta$Group[meta$Group == "Group2"] = "Cell type 2"
  meta$type[meta$type == "cond_A"] = "Condition A"
  meta$type[meta$type == "cond_B"] = "Condition B"
  
  meta = meta[sample(nrow(meta)) , ]
  #title = names(sces)[id]
  #title = str_remove(title , "__n_not_markers_0")
  p = ggplot(meta , aes(x = V1 , y = V2 , col = factor(Group))) + 
    geom_point(alpha = .5 , size = .5) +
    theme_bw() +
    facet_wrap(~type, nrow = 2) + 
    labs(x = "UMAP-1" , y = "UMAP-2") + 
    scale_color_manual(values = cols , name = "Cell type") 
  p
  return(p)
}

3.2 Base simulation

id = 1
p = get_plot(id, "pca.hvgs")
p

ggsave(filename = paste0(figures.dir, "umaps/hvgs/" , names(sces)[id], ".png"), plot = p, width = 4, height = 5)


p = get_plot(id, "pca.all")
p

ggsave(filename = paste0(figures.dir, "umaps/all/" , names(sces)[id], ".png"), plot = p, width = 4, height = 5)

3.3 Only marker perturbed

id = 25
p = get_plot(id, "pca.hvgs")
p

ggsave(filename = paste0(figures.dir, "umaps/hvgs/" , names(sces)[id], ".png"), plot = p, width = 4, height = 5)


p = get_plot(id, "pca.all")
p

ggsave(filename = paste0(figures.dir, "umaps/all/" , names(sces)[id], ".png"), plot = p, width = 4, height = 5)

3.4 Only not markers perturbed

id = 361
p = get_plot(id, "pca.hvgs")
p

ggsave(filename = paste0(figures.dir, "umaps/hvgs/" , names(sces)[id], ".png"), plot = p, width = 4, height = 5)


p = get_plot(id, "pca.all")
p

ggsave(filename = paste0(figures.dir, "umaps/all/" , names(sces)[id], ".png"), plot = p, width = 4, height = 5)

3.5 Both perturbed

id = 385
p = get_plot(id, "pca.hvgs")
p

ggsave(filename = paste0(figures.dir, "umaps/hvgs/" , names(sces)[id], ".png"), plot = p, width = 4, height = 5)


p = get_plot(id, "pca.all")
p

ggsave(filename = paste0(figures.dir, "umaps/all/" , names(sces)[id], ".png"), plot = p, width = 4, height = 5)

3.6 Overall grid (split by sims)

get_all_plots = function(round_sim , red_dim){
  ids = sapply(names(sces) , function(str) str_detect(str , paste0(round_sim, "_n_markers_"), negate = FALSE))
  ids = as.numeric( which(ids) )
  
  plots = bplapply(ids , function(id){
    print(id)
    p = get_plot(id , red_dim)
    return(p)
  }, BPPARAM = mcparam)
  p = ggarrange(plotlist = plots , nrow = 12 , ncol = 12 , common.legend = T)
  return(p)
}


anno = expand.grid(round_sim = c(1,2,3) , red_dim = c("all" , "hvgs"))

plots = lapply(1:nrow(anno) , function(i){
  p = get_all_plots(as.numeric(as.character(anno$round_sim[i])) , red_dim = paste0("pca.", as.character(anno$red_dim[i])))
  ggsave(filename = paste0(figures.dir, "umaps/" , as.character(anno$red_dim[i]), "/" , "all_" , as.character(anno$round_sim[i]), ".png"), plot = p, width = 17, height = 22)
})

4 Session Info

sessionInfo()
## R version 4.2.3 (2023-03-15)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 22.04.2 LTS
## 
## Matrix products: default
## BLAS:   /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3
## LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.20.so
## 
## locale:
##  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
##  [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
##  [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
##  [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
##  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       
## 
## attached base packages:
## [1] stats4    stats     graphics  grDevices utils     datasets  methods  
## [8] base     
## 
## other attached packages:
##  [1] tibble_3.2.1                dplyr_1.1.1                
##  [3] MetBrewer_0.2.0             viridis_0.6.2              
##  [5] viridisLite_0.4.1           ggpubr_0.6.0               
##  [7] ggplot2_3.4.2               stringr_1.5.0              
##  [9] miloDE_0.0.0.9000           miloR_1.6.0                
## [11] edgeR_3.40.2                limma_3.54.2               
## [13] scran_1.26.2                scuttle_1.8.4              
## [15] splatter_1.22.1             geneBasisR_0.0.0.9000      
## [17] BiocParallel_1.32.6         SingleCellExperiment_1.20.1
## [19] SummarizedExperiment_1.28.0 Biobase_2.58.0             
## [21] GenomicRanges_1.50.2        GenomeInfoDb_1.34.9        
## [23] IRanges_2.32.0              S4Vectors_0.36.2           
## [25] BiocGenerics_0.44.0         MatrixGenerics_1.13.0      
## [27] matrixStats_0.63.0          BiocStyle_2.26.0           
## [29] rmarkdown_2.21             
## 
## loaded via a namespace (and not attached):
##   [1] utf8_1.2.3                tidyselect_1.2.0         
##   [3] grid_4.2.3                combinat_0.0-8           
##   [5] munsell_0.5.0             ScaledMatrix_1.6.0       
##   [7] ragg_1.2.5                codetools_0.2-19         
##   [9] statmod_1.5.0             future_1.32.0            
##  [11] tester_0.1.7              withr_2.5.0              
##  [13] batchelor_1.14.1          colorspace_2.1-0         
##  [15] highr_0.10                knitr_1.42               
##  [17] ggsignif_0.6.4            pbmcapply_1.5.1          
##  [19] listenv_0.9.0             labeling_0.4.2           
##  [21] GenomeInfoDbData_1.2.9    mnormt_2.1.1             
##  [23] polyclip_1.10-4           optimParallel_1.0-2      
##  [25] farver_2.1.1              coda_0.19-4              
##  [27] parallelly_1.35.0         vctrs_0.6.2              
##  [29] generics_0.1.3            clusterGeneration_1.3.7  
##  [31] ipred_0.9-14              xfun_0.38                
##  [33] timechange_0.2.0          R6_2.5.1                 
##  [35] doParallel_1.0.17         ggbeeswarm_0.7.1         
##  [37] graphlayouts_0.8.4        rsvd_1.0.5               
##  [39] locfit_1.5-9.7            pals_1.7                 
##  [41] RcppZiggurat_0.1.6        bitops_1.0-7             
##  [43] cachem_1.0.7              DelayedArray_0.24.0      
##  [45] scales_1.2.1              nnet_7.3-18              
##  [47] ggraph_2.1.0              beeswarm_0.4.0           
##  [49] gtable_0.3.3              beachmat_2.14.2          
##  [51] globals_0.16.2            phangorn_2.11.1          
##  [53] tidygraph_1.2.3           timeDate_4022.108        
##  [55] rlang_1.1.0               systemfonts_1.0.4        
##  [57] scatterplot3d_0.3-43      splines_4.2.3            
##  [59] rstatix_0.7.2             yardstick_1.1.0          
##  [61] dichromat_2.0-0.1         broom_1.0.4              
##  [63] checkmate_2.1.0           reshape2_1.4.4           
##  [65] BiocManager_1.30.20       yaml_2.3.7               
##  [67] abind_1.4-5               backports_1.4.1          
##  [69] Rfast_2.0.7               lava_1.7.2.1             
##  [71] tools_4.2.3               bookdown_0.33            
##  [73] Augur_1.0.3               jquerylib_0.1.4          
##  [75] RColorBrewer_1.1-3        plyr_1.8.8               
##  [77] Rcpp_1.0.10               parsnip_1.1.0            
##  [79] sparseMatrixStats_1.11.1  zlibbioc_1.44.0          
##  [81] purrr_1.0.1               RCurl_1.98-1.12          
##  [83] rpart_4.1.19              cowplot_1.1.1            
##  [85] zoo_1.8-12                ggrepel_0.9.3            
##  [87] cluster_2.1.4             furrr_0.3.1              
##  [89] magrittr_2.0.3            magick_2.7.4             
##  [91] data.table_1.14.8         ResidualMatrix_1.8.0     
##  [93] lmtest_0.9-40             patchwork_1.1.2          
##  [95] evaluate_0.20             paleotree_3.4.5          
##  [97] gridExtra_2.3             scater_1.26.1            
##  [99] compiler_4.2.3            maps_3.4.1               
## [101] htmltools_0.5.5           tidyr_1.3.0              
## [103] expm_0.999-7              lubridate_1.9.2          
## [105] DBI_1.1.3                 tweenr_2.0.2             
## [107] MASS_7.3-58.3             Matrix_1.5-4             
## [109] car_3.1-2                 wesanderson_0.3.6        
## [111] cli_3.6.1                 quadprog_1.5-8           
## [113] gdata_2.18.0.1            parallel_4.2.3           
## [115] metapod_1.6.0             gower_1.0.1              
## [117] igraph_1.4.2              pkgconfig_2.0.3          
## [119] RcppGreedySetCover_0.1.0  numDeriv_2016.8-1.1      
## [121] recipes_1.0.5             foreach_1.5.2            
## [123] vipor_0.4.5               bslib_0.4.2              
## [125] ggcorrplot_0.1.4          hardhat_1.3.0            
## [127] dqrng_0.3.0               XVector_0.38.0           
## [129] prodlim_2023.03.31        digest_0.6.31            
## [131] RcppAnnoy_0.0.20          phytools_1.5-1           
## [133] fastmatch_1.1-3           uwot_0.1.14              
## [135] DelayedMatrixStats_1.20.0 gtools_3.9.4             
## [137] lifecycle_1.0.3           nlme_3.1-162             
## [139] jsonlite_1.8.4            carData_3.0-5            
## [141] BiocNeighbors_1.16.0      mapproj_1.2.11           
## [143] fansi_1.0.4               pillar_1.9.0             
## [145] lattice_0.21-8            survival_3.5-5           
## [147] fastmap_1.1.1             plotrix_3.8-2            
## [149] glue_1.6.2                png_0.1-8                
## [151] iterators_1.0.14          bluster_1.8.0            
## [153] ggforce_0.4.1             class_7.3-21             
## [155] stringi_1.7.12            sass_0.4.5               
## [157] textshaping_0.3.6         rsample_1.1.1            
## [159] BiocSingular_1.14.0       future.apply_1.10.0      
## [161] irlba_2.3.5.1             ape_5.7-1